# -*-coding: utf-8 -*-
"""
作者：TianyiFan
日期：2023-年06-月26日
用 途： gurobi 验证
"""

from gurobipy import *
from swo_algo_method import solution_record
import common_fuc
import para as params
import time


class optimal_model:
    def __init__(self,mode = "LP"):
        print("\nGurobi_model 中params.Pair_N = {}".format(params.pair_new))
        self.optimal_rd = solution_record()
        self.solving_time = 0
        self.obj = params.max_bound
        self.solving_time_LB = 0
        self.obj_LB =  params.max_bound
        self.obj_LP =  params.max_bound
        self.model_mode = mode
        self.gurobi_output = []  # 原模型时间构成
    def build_model(self):
        time_start = time.time()
        maxbound =  params.max_bound
        minbound = - params.max_bound
        # region 定义模型
        model = Model("Port_model")
        model.setParam("Timelimit", 300)
        model.setParam("MIPGap", 0.000)
        model.setParam("MIPFocus", 3)
        model.Params.presolve = 1
        # endregion
    
        # region 定义决策变量
        alpha_uih = {}
        alpha_lih = {}
        beta_lik = {}
        beta_uik = {}
        ganma_ik = {}
        for i in range(params.V + 2):
            for h in range(params.H):
                name1 = "alpha_uij_{0:>2}_{1}_".format(i, h)
                name2 = "alpha_lij_{0:>2}_{1}_".format(i, h)
                if self.model_mode =="LP":
                    alpha_uih[i, h] = model.addVar(0, params.c_ui[i], vtype=GRB.CONTINUOUS, name=name1)
                    alpha_lih[i, h] = model.addVar(0, params.c_li[i], vtype=GRB.CONTINUOUS, name=name2)
                else:
                    alpha_uih[i, h] = model.addVar(0, params.c_ui[i], vtype=GRB.INTEGER, name=name1)
                    alpha_lih[i, h] = model.addVar(0, params.c_li[i], vtype=GRB.INTEGER, name=name2)

        for i in range(params.V+ 2):
            for k in range(params.K1):
                name3 = "beta_l{0}_{1}".format(i, k)
                name4 = "beta_u{0}_{1}".format(i, k)
                name5 = "ganma_{0}_{1}".format(i, k)
                if self.model_mode=="LP":
                    beta_lik[i, k] = model.addVar(0, params.z, vtype=GRB.CONTINUOUS, name=name3)
                    beta_uik[i, k] = model.addVar(0, params.z, vtype=GRB.CONTINUOUS, name=name4)
                else:
                    beta_lik[i, k] = model.addVar(0, params.z, vtype=GRB.INTEGER, name=name3)
                    beta_uik[i, k] = model.addVar(0, params.z, vtype=GRB.INTEGER, name=name4)
                ganma_ik[i, k] = model.addVar(0, 1, vtype=GRB.BINARY, name=name5)
        # 辅助决策变量用于线性化
        deta_u_i = {}
        deta_l_i = {}
        abs_ui1 = {}
        abs_li1 = {}
        sum_li1 = {}
        sum_ui1 = {}
        for i in range(1, params.V+ 1):
            name6 = "deta_u_{0}".format(i)
            name7 = "deta_l_{0}".format(i)
            name8 = "abs_u_{0}".format(i)
            name9 = "abs_l_{0}".format(i)
            name10 = "sum_u_{0}".format(i)
            name11 = "sum_i_{0}".format(i)
            if self.model_mode == "LP":
                deta_u_i[i] = model.addVar(0, maxbound, vtype=GRB.CONTINUOUS, name=name6)
                deta_l_i[i] = model.addVar(0, maxbound, vtype=GRB.CONTINUOUS, name=name7)
                abs_ui1[i] = model.addVar(0, maxbound, vtype=GRB.CONTINUOUS, name=name8)
                abs_li1[i] = model.addVar(0, maxbound, vtype=GRB.CONTINUOUS, name=name9)
                sum_ui1[i] = model.addVar(minbound, maxbound, vtype=GRB.CONTINUOUS, name=name10)
                sum_li1[i] = model.addVar(minbound, maxbound, vtype=GRB.CONTINUOUS, name=name11)
            else:
                deta_u_i[i] = model.addVar(0, maxbound, vtype=GRB.INTEGER, name=name6)
                deta_l_i[i] = model.addVar(0, maxbound, vtype=GRB.INTEGER, name=name7)
                abs_ui1[i] = model.addVar(0, maxbound, vtype=GRB.INTEGER, name=name8)
                abs_li1[i] = model.addVar(0, maxbound, vtype=GRB.INTEGER, name=name9)
                sum_ui1[i] = model.addVar(minbound, maxbound, vtype=GRB.INTEGER, name=name10)
                sum_li1[i] = model.addVar(minbound, maxbound, vtype=GRB.INTEGER, name=name11)

        # endregion

        # region 定义目标函数
        obj = LinExpr(0)
        for i in range(1, params.V+ 1):
            expr1 = LinExpr(0)
            expr1.addTerms(params.time_p, deta_l_i[i])
            expr1.addTerms(params.time_p, deta_u_i[i])
            expr1.addTerms(params.time_u, abs_li1[i])
            expr1.addTerms(params.time_u, abs_ui1[i])
            expr2 = LinExpr(0)
            for k in range(params.K1):
                expr2.addTerms(params.t_uki[k][i], beta_uik[i, k])
                expr2.addTerms(params.t_lki[k][i], beta_lik[i, k])
            obj.add(expr1)
            obj.add(expr2)
        # obj.add(quicksum(ganma_ik[i, k] for i in range(params.V + 1) for k in range(params.K1)))
        model.setObjective(obj, GRB.MINIMIZE)
        # endregion

        # region 定义约束条件

        # 目标函数中的等价替换
        for i in range(1, params.V+ 1):
            for h in range(params.H):
                model.addConstr(deta_l_i[i] >= alpha_lih[i, h], name="linearing cons1_{0}_{1}".format(i, h))
                model.addConstr(deta_u_i[i] >= alpha_uih[i, h], name="linearing cons2_{0}_{1}".format(i, h))
        # 最小化目标函数时，绝对值的线性化=恒大于两者最小

        for i in range(1, params.V+ 1):
            expr3 = quicksum(beta_lik[i, k] for k in range(params.K0))
            expr4 = quicksum(beta_uik[i, k] for k in range(params.K0))
            model.addConstr(sum_li1[i] == alpha_lih[i, 0] - expr3, name="Constrain_sum_l_{0}".format(i))
            model.addConstr(sum_ui1[i] == alpha_uih[i, 0] - expr4, name="Constrain_sum_u_{0}".format(i))

        for i in range(1, params.V+ 1):
            model.addGenConstrAbs(abs_li1[i], sum_li1[i], name="abs_constrain_l_{0}".format(i))
            model.addGenConstrAbs(abs_ui1[i], sum_ui1[i], name="abs_constrain_u_{0}".format(i))
        # 约束2/4：每艘船的unloading量守恒、loading量守恒
        if 1 == 1:
            for i in range(1, params.V+ 1):
                expr1 = quicksum(alpha_uih[i,h] for h in range(params.H))
                expr2 = quicksum(alpha_lih[i, h] for h in range(params.H))
                model.addConstr(expr1 == params.c_ui[i], name="Constrain2_u_{0}".format(i))
                model.addConstr(expr2 == params.c_li[i], name="Constrain4_l_{0}".format(i))

        # 约束 3/5：每艘船unloading/loading到堆场的数量守恒
        if 1 == 1:
            for i in range(1, params.V+ 1):
                expr3 = quicksum(beta_uik[i, k] for k in range(params.K1))
                expr4 = quicksum(beta_lik[i, k] for k in range(params.K1))
                model.addConstr(expr3 == params.c_ui[i], name="Constrain_3_u_{0}".format(i))
                model.addConstr(expr4 == params.c_li[i], name="Constrain_5_l_{0}".format(i))

        # 约束6：unloading 和 loading 在每一层的守恒
        if 1 == 1:
            for i in range(1, params.V+ 1):
                for k in range(params.K1):
                    model.addConstr(beta_lik[i, k] + beta_uik[i, k] <= params.z*ganma_ik[i, k], name="Constrain6_{0}_{1}".format(i, k))

        # 约束7-8：每个堆场的loading、unloading 量
        # for i in range(1, params.V+ 1):
        #     for k in range(params.K1):
        #         model.addConstr(beta_lik[i, k] <= params.z * ganma_ik[i, k], name="Constrain7_l_{0}_{1}".format(i, k))
        #         model.addConstr(beta_uik[i, k] <= params.z * ganma_ik[i, k], name="Constrain8_u_{0}_{1}".format(i, k))

        # 约束9：每层的subblock需要小于对应的params.r_ih
        for i in range(1, params.V+ 1):
            expr1 = LinExpr(0)
            expr2 = LinExpr(0)
            for k in range(params.K0):
                expr1.add(ganma_ik[i, k])
            for k in range(params.K0, params.K1):
                expr2.add(ganma_ik[i, k])
            model.addConstr(expr1 <= params.r_ih[i][0], name="Constrain9_0_{0}".format(i))
            model.addConstr(expr2 <= params.r_ih[i][1], name="Constrain9_1_{0}".format(i))

        # 约束10：每个subblock只能分配给一个堆场
        if 1 == 1:
            for k in range(params.K1):
                expr = LinExpr(0)
                for i in range(1, params.V+ 1):
                    expr.add(ganma_ik[i, k])
                model.addConstr(expr <= 1, name="Constrain10_{0}".format(k))

        # 约束11：每个time step，同一层的同一个block内，分配subblock数量不能超过3
        if 1 == 1:
            for t in range(params.T):
                if params.total_handle[t] >= 1:
                    for g in range(params.G):
                        expr = LinExpr(0)
                        for i in range(1,params.V+1):
                        # for i in params.conflict_t[t]:
                            for k in params.block_set[g]:
                                expr.addTerms(params.l_it[i][t], ganma_ik[i, k])
                        model.addConstr(expr <= 3, name="Constrain11_{0}_{1}".format(t, g))

        # 约束12：写法1--从每个subblock看，每个params.neighbor_set中，只能分配给一艘船
        if 1 == 1:
            for t in range(params.T):
                if params.total_handle[t] >= 1:
                    for n in range(params.K1):
                        for k in range(1, len(params.neighbor_set[n])):  # 不包括自己
                            midlist = [n]
                            midlist.append(params.neighbor_set[n][k])
                            expr = LinExpr(0)
                            index_neighbor = n
                            for i in range(1, params.V+ 1):
                            # for i in params.conflict_t[t]:
                                for index in midlist:    # --形成一个pair 对 等价于 for k in neighbor_pair
                                    index_neighbor = index
                                    expr.addTerms(params.l_it[i][t], ganma_ik[i, index])
                            model.addConstr(expr <= 1, name="Constrain12_{0}_({1}-{2})".format(t, n, index_neighbor))

        # 约束13：每个time step，异层转运的数量不能超过一个阈值
        if 1 == 1:
            for t in range(params.T):
                expr = LinExpr(0)
                for i in range(1, params.V+ 1):
                    if params.l_it[i][t] == 1:
                        expr.addTerms(1 / params.handle_i[i], abs_ui1[i])
                        expr.addTerms(1 / params.handle_i[i], abs_li1[i])
                model.addConstr(expr <= params.C, name="Constrain13_{0}".format(t))

        # endregion
        # model.write('double_storey_model.LP')
        model.update()
        model.optimize()
        time_end = time.time()

        print("gurobi_model 中变量数量={}   约束数量={}".format(model.NumVars, model.NumConstrs))
        print("模型的求解状态：", model.getAttr(GRB.attr.Status))
        print('----------optimal status-----------')
        if model.getAttr(GRB.attr.Status) in [2,9]:
            self.solving_time = round(time_end - time_start, 2)
            self.obj = model.ObjVal
            print("模型的最优解 = {1:.2f}  转化为小时={0:.2f} hours".format(model.ObjVal / 3600, model.ObjVal))
            # for i in range(1, params.V+ 1):
            #     print("\nvessel—{} 的求解计划-------------- ".format(i))
            #     print("solution---------输出gamma_k和beta_k的结果-------------")
            #     for k in range(params.K1):
            #         if ganma_ik[i, k].x > 0.5:
            #             print("ganma_{}_{} = 1  ".format(i, k), end=" ")
            #             print("beta_l_{}_{} ={:.2f}  distance={} ".format(i, k, beta_lik[i, k].x, t_lki[k, i]), end=" ")
            #             print("beta_u_{}_{} ={:.2f}  distance={} ".format(i, k, beta_uik[i, k].x, t_uki[k, i]))
            #     print("solution---------------------alpha_uih+alpha_lih")
            #     print("vessel-{0}  params.r_ih=({5}/{6})  上下层unloading: ({1:3.0f},{3:3.0f})  负荷-总负荷：[{7:3.0f}]-[{9}]  loading:({2:3.0f},{4:3.0f})[{8:3.0f}]-[{10}]".
            #           format(i, alpha_uih[i, 0].x, alpha_lih[i, 0].x, alpha_uih[i, 1].x, alpha_lih[i, 1].x,
            #                  params.r_ih[i][0], params.r_ih[i][1], deta_u_i[i].x, deta_l_i[i].x, params.c_ui[i], params.c_li[i]))
            #     print("solution---------------------异层转运数量")
            #     print("vessel-{0}-[{1} ts]  unloading异层转运={2:>3.0f} | {3:>3.0f}|  loading异层转运={4:>3.0f} |{5:>3.0f}| ".
            #           format(i, params.handle_i[i], abs_ui1[i].x, sum_ui1[i].x, abs_li1[i].x, sum_li1[i].x))
            # print("solution---------------------beta_uik+beta_lik")
            # for h in range(params.H):
            #     print("第{0}层的block分配：[堆场号]-(unloading/loading)--------".format(h + 1))
            #     for i in range(1, params.V+ 1):
            #         print("vessel-{0} 的subblock r-ih={1}：".format(i, params.r_ih[i, h]), end=" ")
            #         if h == 0:
            #             for k in range(params.K0):
            #                 if ganma_ik[i, k].x > 0.5:
            #                     print("[{0:>3}]-({1:>3.0f}/{2:>3.0f})".format(k, beta_uik[i, k].x, beta_lik[i, k].x),
            #                           end=" ")
            #             print()
            #         else:
            #             for k in range(params.K0, params.K1):
            #                 if ganma_ik[i, k].x > 0.5:
            #                     print("[{0:>3}]:({1:>3.0f}/{2:>3.0f})".format(k, beta_uik[i, k].x, beta_lik[i, k].x),
            #                           end=" ")
            #             print()
            # 计算每艘船的cost
            cost_t = 0
            cost_p = 0
            cost_u = 0
            cost_i = [0 for _ in range(params.V+ 1)]
            for i in range(1, params.V+ 1):
                cost1 = params.time_p * (max(alpha_uih[i, 0].x, alpha_uih[i, 1].x) + max(alpha_lih[i, 0].x, alpha_lih[i, 1].x))
                unloading_value = 0
                loading_Value = 0
                for k in range(params.K0):
                    unloading_value += beta_uik[i, k].x
                    loading_Value += beta_lik[i, k].x
                # print("abs_ui1 ={}   公式计算的abs_ui1={}".format(abs_ui1[i].x,abs(alpha_uih[i, 0].x - unloading_value)))
                cost2 = params.time_u * (abs(alpha_lih[i, 0].x - loading_Value) + abs(alpha_uih[i, 0].x - unloading_value))
                cost2222 = params.time_u * (abs_ui1[i].x + abs_li1[i].x)
                # print("abs直接计算的结果={}   由原来的公式计算的结果={}".format(cost2222,cost2))
                cost3 = 0
                for k in range(params.K1):
                    cost3 += params.t_lki[k, i] * beta_lik[i, k].x + params.t_uki[k, i] * beta_uik[i, k].x
                cost_i[i] = cost1 + cost2222 + cost3
                cost_t += cost3
                cost_p += cost1
                cost_u += cost2222
                # print("vessel-{} 的总成本={}".format(i, cost_i[i]))
            self.gurobi_output.append(cost_t)
            self.gurobi_output.append(cost_p)
            self.gurobi_output.append(cost_u)
            print("每艘船预留堆场的最优解情况-------------")
            for i in range(1, params.V+ 1):
                print("vessel-{}  cost = {:.2f} 预留subblock = ".format(i,cost_i[i]), end=" ")
                for k in range(params.K1):
                    if ganma_ik[i, k].x > 0.5:
                        print(f"{k:2} ", end="")
                print()
    
            # 将解赋值给solution_record
            gamma_solution = []
            for i in range(1, params.V+ 1):
                self.optimal_rd.cost_i[i] = cost_i[i]
                for h in range(params.H):
                    self.optimal_rd.alpha_lih[i, h] = alpha_lih[i, h].x
                    self.optimal_rd.alpha_uih[i, h] = alpha_uih[i, h].x
                    self.optimal_rd.abs_li1[i] = abs_li1[i].x
                    self.optimal_rd.abs_ui1[i] = abs_ui1[i].x
                temp_solution = []
                for k in range(params.K1):
                    if ganma_ik[i, k].x > 0:
                        temp_solution.append(k)
                        self.optimal_rd.gamma_ik[i, k] = 1
                        self.optimal_rd.beta_lik[i, k] = beta_lik[i, k].x
                        self.optimal_rd.beta_uik[i, k] = beta_uik[i, k].x
                gamma_solution.append(temp_solution)
            params.test_data.data['solution'] = gamma_solution
            # common_fuc.print_pd(params.test_data.data)

        elif model.getAttr(GRB.attr.Status) == 3:
            print("无解")
            print('Optimization was stopped with status =  %d' % model.getAttr(GRB.attr.Status))
            model.computeIIS()
            for c in model.getConstrs():
                if c.IISConstr:
                    print('%s' % c.constrName)
        elif model.getAttr(GRB.attr.Status) == 5:
            print("无界解")

    def build_lower_bound_model(self):
        time_start = time.time()
        maxbound = 999999
        minbound = -999999
        # region 定义模型
        model = Model("Port_model")
        model.setParam("Timelimit", 600)
        model.setParam("MIPGap", 0.000)
        model.setParam("OutputFlag", 0)
        model.Params.presolve = 1
        # endregion

        # region 定义决策变量
        alpha_uih = {}
        alpha_lih = {}
        beta_lik = {}
        beta_uik = {}
        ganma_ik = {}
        for i in range(params.V + 2):
            for h in range(params.H):
                name1 = "alpha_uij_{0:>2}_{1}_".format(i, h)
                name2 = "alpha_lij_{0:>2}_{1}_".format(i, h)
                alpha_uih[i, h] = model.addVar(0, params.c_ui[i], vtype=GRB.CONTINUOUS, name=name1)
                alpha_lih[i, h] = model.addVar(0, params.c_li[i], vtype=GRB.CONTINUOUS, name=name2)

        for i in range(params.V + 2):
            for k in range(params.K1):
                name3 = "beta_l{0}_{1}".format(i, k)
                name4 = "beta_u{0}_{1}".format(i, k)
                name5 = "ganma_{0}_{1}".format(i, k)
                beta_lik[i, k] = model.addVar(0, params.z, vtype=GRB.CONTINUOUS, name=name3)
                beta_uik[i, k] = model.addVar(0, params.z, vtype=GRB.CONTINUOUS, name=name4)
                ganma_ik[i, k] = model.addVar(0, 1, vtype=GRB.CONTINUOUS, name=name5)

        # 辅助决策变量用于线性化
        deta_u_i = {}
        deta_l_i = {}
        abs_ui1 = {}
        abs_li1 = {}
        sum_li1 = {}
        sum_ui1 = {}
        for i in range(1, params.V + 1):
            name6 = "deta_u_{0}".format(i)
            name7 = "deta_l_{0}".format(i)
            name8 = "abs_u_{0}".format(i)
            name9 = "abs_l_{0}".format(i)
            name10 = "sum_u_{0}".format(i)
            name11 = "sum_i_{0}".format(i)
            deta_u_i[i] = model.addVar(0, maxbound, vtype=GRB.CONTINUOUS, name=name6)
            deta_l_i[i] = model.addVar(0, maxbound, vtype=GRB.CONTINUOUS, name=name7)
            abs_ui1[i] = model.addVar(0, maxbound, vtype=GRB.CONTINUOUS, name=name8)
            abs_li1[i] = model.addVar(0, maxbound, vtype=GRB.CONTINUOUS, name=name9)
            sum_ui1[i] = model.addVar(minbound, maxbound, vtype=GRB.CONTINUOUS, name=name10)
            sum_li1[i] = model.addVar(minbound, maxbound, vtype=GRB.CONTINUOUS, name=name11)
        # endregion

        # region 定义目标函数
        obj = LinExpr(0)
        for i in range(1, params.V + 1):
            expr1 = LinExpr(0)
            expr1.addTerms(params.time_p, deta_l_i[i])
            expr1.addTerms(params.time_p, deta_u_i[i])
            expr1.addTerms(params.time_u, abs_li1[i])
            expr1.addTerms(params.time_u, abs_ui1[i])
            expr2 = LinExpr(0)
            for k in range(params.K1):
                expr2.addTerms(params.t_uki[k][i], beta_uik[i, k])
                expr2.addTerms(params.t_lki[k][i], beta_lik[i, k])
            obj.add(expr1)
            obj.add(expr2)
        model.setObjective(obj, GRB.MINIMIZE)
        # endregion

        # region 定义约束条件

        # 目标函数中的等价替换
        for i in range(1, params.V + 1):
            for h in range(params.H):
                model.addConstr(deta_l_i[i] >= alpha_lih[i, h], name="linearing cons1_{0}_{1}".format(i, h))
                model.addConstr(deta_u_i[i] >= alpha_uih[i, h], name="linearing cons2_{0}_{1}".format(i, h))
                # 最小化目标函数时，绝对值的线性化=恒大于两者最小

        for i in range(1, params.V + 1):
            expr3 = quicksum(beta_lik[i, k] for k in range(params.K0))
            expr4 = quicksum(beta_uik[i, k] for k in range(params.K0))
            model.addConstr(sum_li1[i] == alpha_lih[i, 0] - expr3, name="Constrain_sum_l_{0}".format(i))
            model.addConstr(sum_ui1[i] == alpha_uih[i, 0] - expr4, name="Constrain_sum_u_{0}".format(i))

        for i in range(1, params.V + 1):
            model.addGenConstrAbs(abs_li1[i], sum_li1[i], name="abs_constrain_l_{0}".format(i))
            model.addGenConstrAbs(abs_ui1[i], sum_ui1[i], name="abs_constrain_u_{0}".format(i))
        quicksum(beta_lik[i, k] for k in range(params.K0))
        # 约束2/4：每艘船的unloading量守恒、loading量守恒
        if 1 == 1:
            for i in range(1, params.V + 1):
                expr1 = quicksum(alpha_uih[i, h] for h in range(params.H))
                expr2 = quicksum(alpha_lih[i, h] for h in range(params.H))
                model.addConstr(expr1 == params.c_ui[i], name="Constrain2_u_{0}".format(i))
                model.addConstr(expr2 == params.c_li[i], name="Constrain4_l_{0}".format(i))

        # 约束 3/5：每艘船unloading/loading到堆场的数量守恒
        if 1 == 1:
            for i in range(1, params.V + 1):
                expr3 = quicksum(beta_uik[i, k] for k in range(params.K1))
                expr4 = quicksum(beta_lik[i, k] for k in range(params.K1))
                model.addConstr(expr3 == params.c_ui[i], name="Constrain_3_u_{0}".format(i))
                model.addConstr(expr4 == params.c_li[i], name="Constrain_5_l_{0}".format(i))

        # 约束6：unloading 和 loading 在每一层的守恒
        if 1 == 1:
            for i in range(1, params.V + 1):
                for k in range(params.K1):
                    model.addConstr(beta_lik[i, k] + beta_uik[i, k] <= params.z, name="Constrain6_{0}_{1}".format(i, k))

        # 约束7-8：每个堆场的loading、unloading 量
        for i in range(1, params.V + 1):
            for k in range(params.K1):
                model.addConstr(beta_lik[i, k] <= params.z * ganma_ik[i, k], name="Constrain7_l_{0}_{1}".format(i, k))
                model.addConstr(beta_uik[i, k] <= params.z * ganma_ik[i, k], name="Constrain8_u_{0}_{1}".format(i, k))

        # 约束9：每层的subblock需要小于对应的params.r_ih
        for i in range(1, params.V + 1):
            expr1 = LinExpr(0)
            expr2 = LinExpr(0)
            for k in range(params.K0):
                expr1.add(ganma_ik[i, k])
            for k in range(params.K0, params.K1):
                expr2.add(ganma_ik[i, k])
            model.addConstr(expr1 <= params.r_ih[i][0], name="Constrain9_0_{0}".format(i))
            model.addConstr(expr2 <= params.r_ih[i][1], name="Constrain9_1_{0}".format(i))

        # 约束10：每个subblock只能分配给一个堆场
        if 1 == 1:
            for k in range(params.K1):
                expr = LinExpr(0)
                for i in range(1, params.V + 1):
                    expr.add(ganma_ik[i, k])
                model.addConstr(expr <= 1, name="Constrain10_{0}".format(k))

        # 约束11：每个time step，同一层的同一个block内，分配subblock数量不能超过3
        if 1 == 1:
            for t in range(params.T):
                if params.total_handle[t] >= 1:
                    for g in range(params.G):
                        expr = LinExpr(0)
                        for i in range(1, params.V + 1):
                            for k in params.block_set[g]:
                                expr.addTerms(params.l_it[i][t], ganma_ik[i, k])
                        model.addConstr(expr <= 3, name="Constrain11_{0}_{1}".format(t, g))

        # 约束12：写法1--从每个subblock看，每个params.neighbor_set中，只能分配给一艘船
        if 1 == 1:
            for t in range(params.T):
                if params.total_handle[t] >= 1:
                    for n in range(params.K1):
                        for k in range(1, len(params.neighbor_set[n])):  # 不包括自己
                            midlist = [n]
                            midlist.append(params.neighbor_set[n][k])
                            expr = LinExpr(0)
                            index_neighbor = n
                            for i in range(1, params.V + 1):
                                for index in midlist:  # --形成一个pair 对 等价于 for k in neighbor_pair
                                    index_neighbor = index
                                    expr.addTerms(params.l_it[i][t], ganma_ik[i, index])
                            model.addConstr(expr <= 1, name="Constrain12_{0}_({1}-{2})".format(t, n, index_neighbor))

        # 约束13：每个time step，异层转运的数量不能超过一个阈值
        if 1 == 1:
            for t in range(params.T):
                expr = LinExpr(0)
                for i in range(1, params.V + 1):
                    if params.l_it[i][t] == 1:
                        expr.addTerms(1 / params.handle_i[i], abs_ui1[i])
                        expr.addTerms(1 / params.handle_i[i], abs_li1[i])
                model.addConstr(expr <= params.C, name="Constrain13_{0}".format(t))

        # 线性松弛后的补充约束：所有船分配到一个堆场里的数量不超过Z
        for k in range(params.K1):
            expr1 = quicksum(beta_lik[i,k] for i in range(1,params.V+1))
            expr2 = quicksum(beta_uik[i,k] for i in range(1,params.V+1))
            model.addConstr(expr1+expr2 <=params.z,name="constrain14_{}".format(k))

        # endregion
        # veesel_i = [[],[],[],[],[],[]]
        # veesel_i[1] = [23,63]
        # veesel_i[2] = [5,26,28,67,69]
        # veesel_i[3] = [10,15,29,30,32,34,35,72]
        # veesel_i[4] = [38,77]
        # veesel_i[5] = [20,22,24,25,62]
        # for i in range(1,V+1):
        #     for k in veesel_i[i]:
        #         model.addConstr(ganma_ik[i,k] == 1,name="constrain_strict_{}_{}".format(i,k))
        # model.write('double_storey_model.LP')
        model.update()
        model.optimize()
        time_end = time.time()

        print('\n----------optimal status-----------')
        print("gurobi_model 中变量数量={}   约束数量={}".format(model.NumVars, model.NumConstrs))
        print("模型的求解状态：", model.getAttr(GRB.attr.Status))
        if model.getAttr(GRB.attr.Status) in [2, 9]:
            self.solving_time_LB = round(time_end - time_start, 2)
            self.obj_LB = model.ObjVal
            print("Lower_buound 模型的最优解 = {1:.2f}  转化为小时={0:.2f} hours".format(model.ObjVal / 3600, model.ObjVal))
            # print("每艘船预留堆场的最优解情况-------------")
            # for i in range(1, params.V + 1):
            #     for k in range(params.K1):
            #         if ganma_ik[i, k].x > 0:
            #             print("ganma_{}_{} ={:.2f}".format(i, k, ganma_ik[i, k].x))
            # for i in range(1, params.V+ 1):
            #     print("\nvessel—{} 的求解计划-------------- ".format(i))
            #     print("solution---------输出gamma_k和beta_k的结果-------------")
            #     for k in range(params.K1):
            #         if ganma_ik[i, k].x > 0.5:
            #             print("ganma_{}_{} = 1  ".format(i, k), end=" ")
            #             print("beta_l_{}_{} ={:.2f}  distance={} ".format(i, k, beta_lik[i, k].x, t_lki[k, i]), end=" ")
            #             print("beta_u_{}_{} ={:.2f}  distance={} ".format(i, k, beta_uik[i, k].x, t_uki[k, i]))
            #     print("solution---------------------alpha_uih+alpha_lih")
            #     print("vessel-{0}  params.r_ih=({5}/{6})  上下层unloading: ({1:3.0f},{3:3.0f})  负荷-总负荷：[{7:3.0f}]-[{9}]  loading:({2:3.0f},{4:3.0f})[{8:3.0f}]-[{10}]".
            #           format(i, alpha_uih[i, 0].x, alpha_lih[i, 0].x, alpha_uih[i, 1].x, alpha_lih[i, 1].x,
            #                  params.r_ih[i][0], params.r_ih[i][1], deta_u_i[i].x, deta_l_i[i].x, params.c_ui[i], params.c_li[i]))
            #     print("solution---------------------异层转运数量")
            #     print("vessel-{0}-[{1} ts]  unloading异层转运={2:>3.0f} | {3:>3.0f}|  loading异层转运={4:>3.0f} |{5:>3.0f}| ".
            #           format(i, params.handle_i[i], abs_ui1[i].x, sum_ui1[i].x, abs_li1[i].x, sum_li1[i].x))
            # print("solution---------------------beta_uik+beta_lik")
            # for h in range(params.H):
            #     print("第{0}层的block分配：[堆场号]-(unloading/loading)--------".format(h + 1))
            #     for i in range(1, params.V+ 1):
            #         print("vessel-{0} 的subblock r-ih={1}：".format(i, params.r_ih[i, h]), end=" ")
            #         if h == 0:
            #             for k in range(params.K0):
            #                 if ganma_ik[i, k].x > 0.5:
            #                     print("[{0:>3}]-({1:>3.0f}/{2:>3.0f})".format(k, beta_uik[i, k].x, beta_lik[i, k].x),
            #                           end=" ")
            #             print()
            #         else:
            #             for k in range(params.K0, params.K1):
            #                 if ganma_ik[i, k].x > 0.5:
            #                     print("[{0:>3}]:({1:>3.0f}/{2:>3.0f})".format(k, beta_uik[i, k].x, beta_lik[i, k].x),
            #                           end=" ")
            #             print()
            # 计算每艘船的cost
            # cost_i = [0 for _ in range(params.V + 1)]
            # for i in range(1, params.V + 1):
            #     cost1 = params.time_p * (max(alpha_uih[i, 0].x, alpha_uih[i, 1].x) + max(alpha_lih[i, 0].x, alpha_lih[i, 1].x))
            #     unloading_value = 0
            #     loading_Value = 0
            #     for k in range(params.K0):
            #         unloading_value += beta_uik[i, k].x
            #         loading_Value += beta_lik[i, k].x
            #     # print("abs_ui1 ={}   公式计算的abs_ui1={}".format(abs_ui1[i].x,abs(alpha_uih[i, 0].x - unloading_value)))
            #     cost2 = params.time_u * (abs(alpha_lih[i, 0].x - loading_Value) + abs(alpha_uih[i, 0].x - unloading_value))
            #     cost2222 = params.time_u * (abs_ui1[i].x + abs_li1[i].x)
            #     # print("abs直接计算的结果={}   由原来的公式计算的结果={}".format(cost2222,cost2))
            #     cost3 = 0
            #     for k in range(params.K1):
            #         cost3 += params.t_lki[k, i] * beta_lik[i, k].x + params.t_uki[k, i] * beta_uik[i, k].x
            #     cost_i[i] = cost1 + cost2222 + cost3
                # print("vessel-{} 的总成本={}".format(i, cost_i[i]))

            # for i in range(1, params.V + 1):
            #     print("vessel-{}  cost = {:.2f} 预留subblock = ".format(i, cost_i[i]), end=" ")
            #     for k in range(params.K1):
            #         if ganma_ik[i, k].x > 0:
            #             print(f"{k:2} ", end="")
            #     print()
            # 将解赋值给para中记录全局解的静态变量
            # for i in range(1, params.V + 1):
            #     self.optimal_rd.cost_i[i] = cost_i[i]
            #     for h in range(params.H):
            #         self.optimal_rd.alpha_lih[i, h] = alpha_lih[i, h].x
            #         self.optimal_rd.alpha_uih[i, h] = alpha_uih[i, h].x
            #         self.optimal_rd.abs_li1[i] = abs_li1[i].x
            #         self.optimal_rd.abs_ui1[i] = abs_ui1[i].x
            #     for k in range(params.K1):
            #         self.optimal_rd.gamma_ik[i, k] = ganma_ik[i, k].x
            #         self.optimal_rd.beta_lik[i, k] = beta_lik[i, k].x
            #         self.optimal_rd.beta_uik[i, k] = beta_uik[i, k].x

        elif model.getAttr(GRB.attr.Status) == 3:
            print("无解")
            print('Optimization was stopped with status =  %d' % model.getAttr(GRB.attr.Status))
            model.computeIIS()
            for c in model.getConstrs():
                if c.IISConstr:
                    print('%s' % c.constrName)
        elif model.getAttr(GRB.attr.Status) == 5:
            print("无界解")

    def build_relax_model(self):
        time_start = time.time()
        maxbound = 999999
        # region 定义模型
        model = Model("Port_model_Linearing")
        model.setParam("Timelimit", 120)
        model.setParam("MIPGap", 0.000)
        model.setParam("MIPFocus", 3)
        model.setParam("OutputFlag", 0)
        model.Params.presolve = 1
        # endregion

        # region 定义决策变量
        alpha_uih = {}
        alpha_lih = {}
        beta_lik = {}
        beta_uik = {}
        ganma_ik = {}
        for i in range(params.V + 2):
            for h in range(params.H):
                name1 = "alpha_uij_{0:>2}_{1}_".format(i, h)
                name2 = "alpha_lij_{0:>2}_{1}_".format(i, h)
                alpha_uih[i, h] = model.addVar(0, params.c_ui[i], vtype=GRB.INTEGER, name=name1)
                alpha_lih[i, h] = model.addVar(0, params.c_li[i], vtype=GRB.INTEGER, name=name2)

        for i in range(params.V + 2):
            for k in range(params.K1):
                name3 = "beta_l{0}_{1}".format(i, k)
                name4 = "beta_u{0}_{1}".format(i, k)
                name5 = "ganma_{0}_{1}".format(i, k)
                beta_lik[i, k] = model.addVar(0, params.z, vtype=GRB.INTEGER, name=name3)
                beta_uik[i, k] = model.addVar(0, params.z, vtype=GRB.INTEGER, name=name4)
                ganma_ik[i, k] = model.addVar(0, 1, vtype=GRB.BINARY, name=name5)

        # 辅助决策变量用于线性化
        deta_u_i = {}
        deta_l_i = {}
        abs_ui1 = {}
        abs_li1 = {}
        for i in range(1, params.V + 1):
            name6 = "deta_u_{0}".format(i)
            name7 = "deta_l_{0}".format(i)
            name8 = "abs_u_{0}".format(i)
            name9 = "abs_l_{0}".format(i)
            deta_u_i[i] = model.addVar(0, maxbound, vtype=GRB.INTEGER, name=name6)
            deta_l_i[i] = model.addVar(0, maxbound, vtype=GRB.INTEGER, name=name7)
            abs_ui1[i] = model.addVar(0, maxbound, vtype=GRB.INTEGER, name=name8)
            abs_li1[i] = model.addVar(0, maxbound, vtype=GRB.INTEGER, name=name9)
        # endregion

        # region 定义目标函数
        obj = LinExpr(0)
        for i in range(1, params.V + 1):
            expr1 = LinExpr(0)
            expr1.addTerms(params.time_p, deta_l_i[i])
            expr1.addTerms(params.time_p, deta_u_i[i])
            expr1.addTerms(params.time_u, abs_li1[i])
            expr1.addTerms(params.time_u, abs_ui1[i])
            expr2 = LinExpr(0)
            for k in range(params.K1):
                expr2.addTerms(params.t_uki[k][i], beta_uik[i, k])
                expr2.addTerms(params.t_lki[k][i], beta_lik[i, k])
            obj.add(expr1)
            obj.add(expr2)
        model.setObjective(obj, GRB.MINIMIZE)
        # endregion

        # region 定义约束条件

        # 目标函数中的等价替换
        for i in range(1, params.V + 1):
            for h in range(params.H):
                model.addConstr(deta_l_i[i] >= alpha_lih[i, h], name="linearing cons1_{0}_{1}".format(i, h))
                model.addConstr(deta_u_i[i] >= alpha_uih[i, h], name="linearing cons2_{0}_{1}".format(i, h))
        # 最小化目标函数时，绝对值的线性化=恒大于两者最小

        for i in range(1, params.V + 1):
            expr3 = quicksum(beta_lik[i, k] for k in range(params.K0))
            expr4 = quicksum(beta_uik[i, k] for k in range(params.K0))
            model.addConstr(abs_li1[i] >= alpha_lih[i, 0] - expr3, name="Constrain_absl_plus_{0}".format(i))
            model.addConstr(abs_li1[i] >= expr3 - alpha_lih[i, 0], name="Constrain_absl_diff_{0}".format(i))
            model.addConstr(abs_ui1[i] >= alpha_uih[i, 0] - expr4, name="Constrain_absu_plus_{0}".format(i))
            model.addConstr(abs_ui1[i] >= expr4 - alpha_uih[i, 0], name="Constrain_absu_diff_{0}".format(i))

        # 约束2/4：每艘船的unloading量守恒、loading量守恒
        if 1 == 1:
            for i in range(1, params.V + 1):
                expr1 = quicksum(alpha_uih[i, h] for h in range(params.H))
                expr2 = quicksum(alpha_lih[i, h] for h in range(params.H))
                model.addConstr(expr1 == params.c_ui[i], name="Constrain2_u_{0}".format(i))
                model.addConstr(expr2 == params.c_li[i], name="Constrain4_l_{0}".format(i))

        # 约束 3/5：每艘船unloading/loading到堆场的数量守恒
        if 1 == 1:
            for i in range(1, params.V + 1):
                expr3 = quicksum(beta_uik[i, k] for k in range(params.K1))
                expr4 = quicksum(beta_lik[i, k] for k in range(params.K1))
                model.addConstr(expr3 == params.c_ui[i], name="Constrain_3_u_{0}".format(i))
                model.addConstr(expr4 == params.c_li[i], name="Constrain_5_l_{0}".format(i))

        # 约束6：unloading 和 loading 在每一层的守恒
        if 1 == 1:
            for i in range(1, params.V + 1):
                for k in range(params.K1):
                    model.addConstr(beta_lik[i, k] + beta_uik[i, k] <= params.z, name="Constrain6_{0}_{1}".format(i, k))

        # 约束7-8：每个堆场的loading、unloading 量
        for i in range(1, params.V + 1):
            for k in range(params.K1):
                model.addConstr(beta_lik[i, k] <= params.z * ganma_ik[i, k], name="Constrain7_l_{0}_{1}".format(i, k))
                model.addConstr(beta_uik[i, k] <= params.z * ganma_ik[i, k], name="Constrain8_u_{0}_{1}".format(i, k))

        # 约束9：每层的subblock需要小于对应的params.r_ih
        for i in range(1, params.V + 1):
            expr1 = LinExpr(0)
            expr2 = LinExpr(0)
            for k in range(params.K0):
                expr1.add(ganma_ik[i, k])
            for k in range(params.K0, params.K1):
                expr2.add(ganma_ik[i, k])
            model.addConstr(expr1 <= params.r_ih[i][0], name="Constrain9_0_{0}".format(i))
            model.addConstr(expr2 <= params.r_ih[i][1], name="Constrain9_1_{0}".format(i))

        # 约束10：每个subblock只能分配给一个堆场
        if 1 == 1:
            for k in range(params.K1):
                expr = LinExpr(0)
                for i in range(1, params.V + 1):
                    expr.add(ganma_ik[i, k])
                model.addConstr(expr <= 1, name="Constrain10_{0}".format(k))

        # 约束11：每个time step，同一层的同一个block内，分配subblock数量不能超过3
        if 1 == 1:
            for t in range(params.T):
                if params.total_handle[t] >= 1:
                    for g in range(params.G):
                        expr = LinExpr(0)
                        for i in range(1, params.V + 1):
                            for k in params.block_set[g]:
                                expr.addTerms(params.l_it[i][t], ganma_ik[i, k])
                        model.addConstr(expr <= 3, name="Constrain11_{0}_{1}".format(t, g))

        # 约束12：写法1--从每个subblock看，每个params.neighbor_set中，只能分配给一艘船
        if 1 == 1:
            for t in range(params.T):
                if params.total_handle[t] >= 1:
                    for n in range(params.K1):
                        for k in range(1, len(params.neighbor_set[n])):  # 不包括自己
                            midlist = [n]
                            midlist.append(params.neighbor_set[n][k])
                            expr = LinExpr(0)
                            index_neighbor = n
                            for i in range(1, params.V + 1):
                                for index in midlist:  # --形成一个pair 对 等价于 for k in neighbor_pair
                                    index_neighbor = index
                                    expr.addTerms(params.l_it[i][t], ganma_ik[i, index])
                            model.addConstr(expr <= 1, name="Constrain12_{0}_({1}-{2})".format(t, n, index_neighbor))

        # 约束13：每个time step，异层转运的数量不能超过一个阈值
        if 1 == 1:
            for t in range(params.T):
                expr = LinExpr(0)
                for i in range(1, params.V + 1):
                    if params.l_it[i][t] == 1:
                        expr.addTerms(1 / params.handle_i[i], abs_ui1[i])
                        expr.addTerms(1 / params.handle_i[i], abs_li1[i])
                model.addConstr(expr <= params.C, name="Constrain13_{0}".format(t))

        # endregion
        # veesel_i = [[],[],[],[],[],[]]
        # veesel_i[1] = [23,63]
        # veesel_i[2] = [5,26,28,67,69]
        # veesel_i[3] = [10,15,29,30,32,34,35,72]
        # veesel_i[4] = [38,77]
        # veesel_i[5] = [20,22,24,25,62]
        # for i in range(1,V+1):
        #     for k in veesel_i[i]:
        #         model.addConstr(ganma_ik[i,k] == 1,name="constrain_strict_{}_{}".format(i,k))
        # model.write('double_storey_model.LP')
        model.update()
        model.optimize()
        time_end = time.time()

        print("gurobi_model 中变量数量={}   约束数量={}".format(model.NumVars, model.NumConstrs))
        print("模型的求解状态：", model.getAttr(GRB.attr.Status))
        print('\n----------optimal status-----------')
        if model.status != GRB.INFEASIBLE and model.status != GRB.UNBOUNDED:
            self.solving_time = round(time_end - time_start, 2)
            self.obj_LP = model.ObjVal
            print("LP模型的最优解 = {1:.2f}  转化为小时={0:.2f} hours".format(model.ObjVal / 3600, model.ObjVal))

            # 计算每艘船的cost
            cost_i = [0 for _ in range(params.V + 1)]
            for i in range(1, params.V + 1):
                cost1 = params.time_p * (max(alpha_uih[i, 0].x, alpha_uih[i, 1].x) + max(alpha_lih[i, 0].x, alpha_lih[i, 1].x))
                unloading_value = 0
                loading_Value = 0
                for k in range(params.K0):
                    unloading_value += beta_uik[i, k].x
                    loading_Value += beta_lik[i, k].x
                # print("abs_ui1 ={}   公式计算的abs_ui1={}".format(abs_ui1[i].x,abs(alpha_uih[i, 0].x - unloading_value)))
                cost2 = params.time_u * (abs(alpha_lih[i, 0].x - loading_Value) + abs(alpha_uih[i, 0].x - unloading_value))
                cost2222 = params.time_u * (abs_ui1[i].x + abs_li1[i].x)
                # print("abs直接计算的结果={}   由原来的公式计算的结果={}".format(cost2222,cost2))
                cost3 = 0
                for k in range(params.K1):
                    cost3 += params.t_lki[k, i] * beta_lik[i, k].x + params.t_uki[k, i] * beta_uik[i, k].x
                cost_i[i] = cost1 + cost2222 + cost3
                # print("vessel-{} 的总成本={}".format(i, cost_i[i]))
            print("每艘船预留堆场的最优解情况-------------")
            for i in range(1, params.V + 1):
                print("vessel-{}  cost = {:.2f} 预留subblock = ".format(i, cost_i[i]), end=" ")
                for k in range(params.K1):
                    if ganma_ik[i, k].x > 0.5:
                        print(f"{k:2} ", end="")
                print()

            # # 将解赋值给para中记录全局解的静态变量
            # for i in range(1, params.V + 1):
            #     self.optimal_rd.cost_i[i] = cost_i[i]
            #     for h in range(params.H):
            #         self.optimal_rd.alpha_lih[i, h] = alpha_lih[i, h].x
            #         self.optimal_rd.alpha_uih[i, h] = alpha_uih[i, h].x
            #         self.optimal_rd.abs_li1[i] = abs_li1[i].x
            #         self.optimal_rd.abs_ui1[i] = abs_ui1[i].x
            #     for k in range(params.K1):
            #         self.optimal_rd.gamma_ik[i, k] = ganma_ik[i, k].x
            #         self.optimal_rd.beta_lik[i, k] = beta_lik[i, k].x
            #         self.optimal_rd.beta_uik[i, k] = beta_uik[i, k].x
        elif model.getAttr(GRB.attr.Status) == 3:
            print("无解")
            print('Optimization was stopped with status =  %d' % model.getAttr(GRB.attr.Status))
            model.computeIIS()
            for c in model.getConstrs():
                if c.IISConstr:
                    print('%s' % c.constrName)
        elif model.getAttr(GRB.attr.Status) == 5:
            print("无界解")

